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Abstract 

A numerical simulation is performed of the gravitational collapse of a spher- 
ically symmetric scalar field. The algorithm uses the null initial value for- 
mulation of the Einstein-scalar equations, but does not use adaptive mesh 
refinement. A study is made of the critical phenomena found by Choptuik 
in this system. In particular it is verified that the critical solution exhibits 
periodic self-similarity. This work thus provides a simple algorithm that gives 
verification of the Choptuik results. 
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I. INTRODUCTION 



Recently Choptuik has discovered scaling behavior in the collapse of a spherically sym- 
metric scalar field to form a black hole. Choptuik numerically evolves a family of initial 
data parametrized by p to find that the mass of the black hole is M oc — p*y where p* 
is the critical value of p and 7 is the scaling exponent. For the data with p = p*, a zero 
mass singularity forms. This critical solution has the property of periodic self-similarity: 
the scalar field evolves, after a certain amount of time, to a copy of its profile with the scale 
of space shrunk. Similar results have been found by Abrahams and Evans for vacuum 
axisymmetric gravitational collapse. 

To treat the critical solution the parameter p must be tuned to p* to great accuracy. 
In addition the size of features must be resolved on extremely small scales. Therefore one 
might worry that the periodic self-similarity of the critical solution could be an artifact of the 
numerical algorithms used rather than the actual behavior of the collapse of a scalar field. 
The most srtaightforward way to show that the results of [jl| are not numerical artifacts is 
to perform a numerical treatment of the same physical problem using a completely different 
algorithm. However, any accurate treatment of the critical solution must be able to resolve 
features on extremely small scales. Choptuik achieved this by using an adaptive mesh 
refinement algorithm. Since adaptive mesh refinement algorithms are fairly complicated it 
is not a trivial task to produce another adaptive mesh refinement code to redo the Choptuik 
result. Even with such a code one might worry that the results are an artifact of adaptive 
mesh refinement. In this paper I present results of a numerical simulation of the critical 
collapse of a scalar field. The algorithm uses the null initial value formulation of the problem, 
rather than the spacelike initial value formulation used in . In addition, no adaptive mesh 
refinement is used. I find results in agreement with those of Choptuik. Section II describes 
the null initial formulation of the collapse of a spherically symmetric scalar field. Section III 
is a description of the numerical algorithm used to evolve this system. Section IV contains 
the results. 
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II. NULL INITIAL VALUE FORMULATION 



The null initial value formulation for the collapse of a spherically symmetric scalar field 
was worked out by Christodoulou. In this section we review this formulation and introduce 
the notation to be used in the rest of the paper. The Einstein-scalar equations are 

Rab = 87rV,$V,$ (1) 

where Rab is the Ricci tensor and $ is the scalar field. 

Null coordinates u and v are defined as follows: u is the proper time of an observer at 
the origin and is constant on outgoing light rays; m = is the initial data surface. The 
coordinate v is constant on ingoing light rays and is equal to the usual area coordinate r on 
the initial data surface. (In what follows r will be regarded as a function of u and v). 

For any quantity / define /, /' and / by 

. (3) 

-If" 

f=- f{u,v)r'{u,v)dv . (4) 



r JO 



Let h be the scalar such that /i = $. Then as a consequence of Einstein's equations the 
metric can be written in terms of h and r. Define the quantities q and g by 

q^r-'{h - hf , (5) 

g = exp (Anrq) . (6) 

Then the metric is 

ds'^ = -2gr' dudv + P dn'^ . (7) 
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Here dfl^ is the two-sphere metric. Einstein's equations also provide evolution equations for 
h and r These are HI 



2r 



l:i9--9){h-h) , (8) 



r = - ^9 ■ (9) 



III. NUMERICAL METHODS 

A numerical simulation of these equations was first performed by Goldwirth and Piran. 
Q A version of this algorithm was applied to the Choptuik problem by Grundlach, Price 
and PuUin. |Q However, the methods of references are not accurate enough to treat the 
critical solution. I will start by discussing those features that my algorithm has in common 
with these earlier treatments and then present the new features that give the improvements 
in accuracy. 

Initial data for the Einstein-scalar equations is just the value of h on the initial data sur- 
face. One evolves these equations as follows: first one finds in succession the quantities 
h, q, 9 and g. This is done by evaluating the integrals for these quantities using Simpson's 
rule for unequally spaced points. Then the equations for h and f are used to evolve h and 
r forward one time step. This process is iterated as many times as necessary: i.e. until 
either a black hole forms or the field disperses. To find whether a black hole forms one looks 
for a marginally outer trapped surface. For spherical symmetry this is a surface for which 
V^rVa?" = 0. (In practice the code cannot evolve up to the marginally outer trapped surface 
so one looks for the condition V^rVar ^0.) Since such a surface has r = 2M the mass of 
the black hole is then half the radius of the marginally outer trapped surface. 

In this method was used to study the scaling behavior of the mass of the black hole. 
However, the method is not accurate enough for a treatment of the critical solution. I will 
now describe the sources of inaccuracy and the methods that my code uses to overcome 
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them. One source of inaccuracy comes from the fact that one divides by r in evaluating the 
quantities h, q and g. This leads to inaccuracies near r = 0. The second source of inaccuracy 
comes from the behavior of the critical solution. As the critical solution evolves its structure 
appears on ever smaller spatial scales. With a fixed spatial resolution there comes a time 
when the number of grid points is not sufficient to resolve the structure of the scalar field. 

My code overcomes the first source of inaccuracy as follows: first expand h in a Taylor 
series in r. 

h^ho + hir + O(r') , (10) 
Then the Taylor series for h, q, g and g are given by 

h^ho + ^hr + O(r^) , (11) 

q^\hlr + 0{r') , (12) 
g^l + ^hlr + O(r^) . (13) 

g^l + Ihjr + 0{r^) . (14) 

Thus to find the behavior of all quantities near the origin one needs to find only hg and hi. 
This is done by fitting the first four values of /i to a hne. Equations (11-14) are then used to 
find the values of h, q, g and g for the first three values of r. The Simpson's rule integration 
is then used for all other values of r. 

The inaccuracy due to poor spatial resolution is solved as follows: first note that the 
inaccuracy is due to the size of the spatial structure being much smaller than the overall 
size of the grid. However, in the null initial value formulation the grid points are tied to 
ingoing hght rays. Thus as the system evolves the overall size of the grid becomes smaller 
and the resolution improves. The critical solution forms a zero mass singularity. Consider 
the ingoing light ray that just barely hits this singularity. Choose the outermost gridpoint 
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to correspond to this light ray. Then the size of the grid shrinks as much as the size of 
the spatial features and thus the resolution is always good. In practice one does not know 
beforehand the position on the initial data surface of the light ray that just barely hits the 
singularity. Therefore I have run the code as follows: First I make an estimate of where the 
outermost gridpoint should be and choose the point to be slightly higher than the estimate. 
This ensures that the evolution will proceed for a while with good resoltion. However, after 
a certain amount of time, because the outermost grid point is too far, the spatial features 
will come to occupy a relatively small number of grid points. Observation of the size of 
the spatial features at this time allows a refined estimate for the position of the outermost 
grid point. This refined estimate is then used to choose a (shghtly too large) outermost 
gridpoint. This whole process is then iterated (a few times) until the outermost gridpoint 
corresponds, with good accuracy to the light ray that just barely hits the singularity. 

In this way the overall size of the grid shrinks in proportion to the size of the spatial 
features. However, it is still the case that the number of grid points decreases during the 
evolution. This is because each grid point corresponds to an ingoing light ray. The grid 
point is therefore lost when the hght ray hits the origin. This difficulty is dealt with in the 
code as follows: the evolution proceeds until half of the grid points are lost. These grid 
points are then interpolated halfway in between each of the remaining grid points. Thus 
over time the number of gridpoints is maintained. 

This algorithm is far simpler (though far more specialized) than adaptive mesh refine- 
ment. The code was about 200 lines of Fortran. 

Although this algorithm is well suited for a study of the critical solution, it is not so well 
suited for a study of the scaling behavior of the mass of the black hole. In treating the critical 
solution a great deal of accuracy was gained by choosing the outermost gridpoint to be the 
null geodesic that just barely hits the singularity. However, in treating the formation of a 
black hole the outermost gridpoint must be chosen sufficiently far that the corresponding 
geodesic does not hit the origin before the formation of a marginally outer trapped surface. 
When treating a range of p the outermost gridpoint must satisfy this condition for every 
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spacetime in the one parameter family. But then for p sufficiently close to p* the spatial 
resolution will not be good enough to resolve the spatial structure of the solution. Therefore 
only a limitted range of p can be treated using this algorithm. 



All runs were done with 300 spatial gridpoints. The code was developed and run on a 
NeXT workstation. Then for higher accuracy it was run on a Cray Y-MP. The initial data 
for the scalar field $ on the u — Q surface was 



where 0o, ''o and a are constants. Thus the initial value of $ has essentially a gaussian 
profile of width a centered about a spherical shell of radius Tq. The values of a and Tq are 
fixed, while the amplitude 0o is the parameter p. 

The code was first run to find the critical value of the parameter p. This was done by 
a binary search: a parameter that led to the formation of a black hole was higher than 
the critical parameter. One that did not lead to black hole formation was smaller than the 
critical value. The next estimate of the critical parameter was chosen halfway between one 
known to be too high and one known to be too low. This process was iterated until the 
critical value of the parameter was found to the needed accuracy. At the same time the 
value of the outermost gridpoint was chosen as described in the previous section. 

The code was then run to examine the scaling behavior of the black hole mass. The 
outermost gridpoint was chosen large enough to accomodate a range of the parameter p to 
good accuracy. The code was then run for several values of p in this range. For each value 
of p the code was run until a marginally outer trapped surface formed. The mass M of the 
black hole was then calculated. In figure 1 In M is plotted as a function of In |p — p* | . 



IV. RESULTS 




(15) 
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Figure 1 

The points are well fit by a straight line whose slope is ~ 0.38. Thus the mass behaves 
like M (X {p — p*y where 7 ~ 0.38. This is consistent with the results of references |T|J§]. 

Next the code was run to treat the critical solution. The parameter p was set to p* and 
the outermost gridpoint was set to its optimum value. Let u* be the value of u at which the 
singularity forms. Define T and i? by T = — In {u* — u) and R = re^. Then the periodic 
self-similar property of the critical solution is that h{R, T) is a periodic function of T. 
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Figure 2 

In figure 2 the quantity h is plotted as a function of R and T. After a certain amount 
of evolution the scalar field settles down to a behavior that seems to be periodic in T. To 
examine this apparent periodicity more carefully we pick an identifiable set of times: those 
times at which the maximum of h occurs at r = 0. The corresponding values of T are 
Ti 2.58 , T2 ^ 6.02 , T3 f=:J 9.47 and T4 12.95. Note that these times are equally spaced 
in T with a spacing AT = 3.45 in agreement with the result found by Choptuik. We then 
plot h{R) at each of these times to see whether this function is the same at each time. 
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Figure 3 

Figure 3 shows h{R,Ti). 
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Figure 4 

Figure 4 shows h{R,Ti) and T2) Note that the functions agree. I emphasize that 
this figure has two different functions h{R,Ti) and h{R,T2) plotted on the same graph. 
The agreement is so good that one cannot tell that the two functions differ. 
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Figure 5 

Figure 5 shows h{R, Ti) , h{R, T2) and h{R, T^) plotted together on the same graph. 
Again the agreement is so good that one cannot tell that the functions differ. 
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Figure 6 

Figure 6 shows h{R,Ti) , T2) , h{R,T^) and hi^R^T^) plotted together on the same 
graph. Here the disagreement between the functions is barely visible. 

Thus it is clear that after some initial evolution the scalar field of the critical solution 
settles down to an evolution that is periodic in T. Therefore we have confirmed, using a 
completely different algorithm from that of reference 0], that the critical solution spacetime 
has periodic self-similarity. 
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FIGURES 

FIG. 1. Scaling of the black hole mass. InM is plotted vs. In (p — p*)- The result is a straight 
line whose slope is the critical exponent 7 « 0.38. 

FIG. 2. Behavior of /i as a function of the logarithmic coordinates R and T. After some ini- 
tial evolution /i is a periodic function of T. This demonstrates the periodic self-similarity of the 
spacetime. 

FIG. 3. h{R) is plotted at Ti 

FIG. 4. h{R) is plotted at two different times T\ and T2 on the same graph. Note that the two 
functions agree completely. 

FIG. 5. h{R) is plotted at three different times Ti, T2 and T3 on the same graph. Note that 
the three functions agree completely. 

FIG. 6. h{R) is plotted at four different times Ti, T2, T3 and T4 on the same graph. One can 
barely see the disagreement among these functions. 
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